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Premelting describes the confluence of phenomena that are responsible for the stable existence 
of the liquid phase of matter in the solid region of its bulk phase diagram. Here we develop a 
theoretical description of the premelting of water ice contained in a porous matrix, made of a 
material with a melting temperature substantially larger than ice itself, to predict the amount of 
liquid water in the matrix at temperatures below its bulk freezing point. Our theory combines the 
interfacial premelting of ice in contact with the matrix, grain boundary melting in the ice, and 
impurity and curvature induced premelting, the latter occurring in regions which force the ice-liquid 
interface into a high curvature configuration. These regions are typically found at points where the 
matrix surface is concave, along contact lines of a grain boundary with the matrix, and in liquid 
veins. Both interfacial premelting and curvature induced premelting depend on the concentration 
of impurities in the liquid, which, due to the small segregation coefficient of impurities in ice are 
treated as homogeneously distributed in the premelted liquid. Our principal result is an equation 
for the fraction of liquid in the porous medium as a function of the undercooling, which embodies 
the combined effects of interfacial premelting, curvature induced premelting, and impurities. The 
result is analyzed in detail and applied to a range of experimentally relevant settings. 



I. INTRODUCTION 

Due to advances in our understanding of the pre- 
melting, of ice and other materials over the last fifteen 
years and an increasing appreciation of the influence 
and extent of this phenomenon in technological, terres- 
trial and extraterrestrial settings, it is prudent to revisit 
the detailed thermodynamic basis of its manifestation in 
porous media where water is invariably found in adsorbed 
layers or in bulk. In 1992 Cahn, Dash, and Fu (CDF) 
[2| developed the original premelting based theory for 
the liquid fraction /; of water in a porous medium at 
sub-freezing temperatures. Given experimental and the- 
oretical advances in the interim, we extend their work 
in a number of respects. CDF considered a matrix of 
monodisperse spheres arranged either in a simple cubic 
packing or in a cubic close packing, regarded as lower and 
upper bounds for the experimental systems. We design 
our theory for the random close packing of monodisperse 
spheres, which may be a more realistic description for 
most applications, e.g. frozen soils, where even a small 
degree of polydispersity prevents a periodic packing of 
the matrix particles. The assumed periodic packing of 
the matrix particles leads CDF to take the grain size 
of the interstitial ice to be on the order of the size of 
the matrix lattice unit cell. Given that the degree of 
polycrystallinity of the interstitial ice has a complicated 
dependence on the freezing rate [see e.g., Ref. [H, and ref- 
erences therein] and the sample annealing history, the ice 
grain size may differ substantially from the matrix parti- 
cle size. Whence, polycrystallinity of the interstitial ice 
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enters into our theory as a parameter, namely the den- 
sity of grain boundary surface relative to the interstitial 
volume. Calculations since CDF showed that the dis- 
persion force contribution to interfacial melting could be 
repulsive or attractive [l], Q and thus favor or suppress 
melting in a manner and strength that depends on the 
effective dielectric response of the ice/water/matrix in- 
terface. Moreover, it has been understood that the great 
inter-experimental variability in both the interfacial and 
surface melting of ice was likely due to the sensitivity 
of the phenomena to impurities which have a rich and 
complex influence on the melting behavior p|. Thus, our 
theory allows for both attractive and repulsive dispersion 
forces of arbitrary strength (depending on the matrix ma- 
terial) and includes the effect of charges at the interface, 
which give rise to an electrostatic force. Finally, most 
importantly, our work includes the effect of charged im- 
purities in the premelted liquid on both interfacial and 
curvature induced premelting. 



The article is structured as follows. In Section [H] we 
lay out the tenets and details of our theory, discussing in- 
terfacial premelting and curvature induced premelting in 
the presence of charged impurities and combine them in 
order to describe premelting in a porous medium. The fi- 
nal result of the theoretical treatment is an implicit func- 
tional equation for the liquid fraction fi, Eq. (|2"Tj) . The 
first part of Section Hill is dedicated to the analysis of in- 
terfacial premelting as described by Eq. (|7|), and in the 
second part Eq. (f2~Tj) for the liquid fraction // is analyzed 
by extracting exact solutions for certain limiting cases 
and numerical solution in the general case. The porous 
media which we consider consist of fused quartz, gold, 
and silicon. We summarize in Section [IVI 
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II. THEORY 



Interfacial melting and grain boundary melting: 
Planar interfaces 



In order to calculate the equilibrium thickness d of a 
premelted layer for a given undercooling AT = T m — T, 
where T m is the bulk melting temperature of ice and T 
is the temperature of the system, we consider a three 
layer planar system. Layer / is given by the substrate, 
i.e., in the case of interfacial melting, the layer consists 
of the material which forms the porous medium, or of ice 
when grain boundary melting is considered. Layer II is 
the premelted quasi-liquid layer of thickness d and finally 
layer III consists of ice. In what follows we consider the 
complete free energy of the system and then focus on 
the individual complications and contributions in turn 
to show how known asymptotic results emerge from the 
theory. 



1. General theory in the presence of impurities 

It is assumed that layer II contains a certain amount 
of impurities. The Gibbs free energy of the system per 
unit area is given by 

G(T,P,N 3 ,N h Ni) 

= ii a (T,P)N a +in(T,P)Ni+m(T,P)Ni 



-R g TN l 



In — - 1 



-2(d) 



(1) 



where T and P are temperature and pressure, p, s , fii, 
Hi are the chemical potentials per mole of the solid, the 
liquid, and the impurity species, respectively, N s , Ni, 
Ni are the numbers of mole per unit area of the solid, 
the liquid, and the impurities, respectively. The term 

RgTNi 
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where R g is the gas constant, corre- 
sponds to the mixing entropy of the impurities assuming 
low impurity concentration so that the ideal gas result 
applies. Finally, 1(d) denotes the total effective interfa- 
cial free energy described in more detail below. Note that 
d is directly proportional to fVj through Ni = pid, where 
pi is the molar density of the liquid which is a constant 
in the calculation. 

The solid (ice) and the liquid (water) parts of the sys- 
tem can exchange material and hence the condition for 
thermodynamic equilibrium is given by 

0G_ = dG_ 

dN s dNi' 1 ' 

Applying this criterion to the free energy given above, 
and using the fact that dp St i/dT = — s Si ; where s s (s;) is 
the specific entropy of the solid (liquid), we can approx- 
imate pi — fi s ~ q m Tjr-, where q m is the latent heat of 
melting per mole and arrive at 

AT R g T m Nj 
PiQmT^- = — —, T(d). (3) 



In the absence of interfacial forces we have I' (d) = 0, and 
the result reduces to AT = q s ™ pi, where the molar 
density of impurities is pi = Ni/d. This undercooling 
corresponds precisely to the classical colligativc effect- 
the depression of the bulk freezing temperature due to 
dissolved impurities. Obviously, compared to the force 
free case, repulsive interfacial forces for which I'(d) < 
result in a larger undercooling AT thereby enhancing the 
layer thickness. Conversely, attractive interfacial forces 
wherein I'(d) > suppress the formation of a premelted 
liquid film. Indeed, for sufficiently large AT, Eq. @ may 
have no solution showing that the attractive interfacial 
forces can dominate the colligativc effect of the impurities 
and the premelted liquid layer will vanish. 



2. Contributions to the effective interfacial free energy 1(d) 

The interfacial energy has two contributions; T(d) = 
Fdis(d) + F e \ ec (d), where Fdi s (d) is due to dispersion 
forces, acting between the three interfaces defining the 
system, and F e \ ec (d) captures the effect of immobile 
charges at the interface. Although we compute the 
strength of the dispersion force interaction from the full 
frequency dependent theory we take the small thick- 
ness limit in which only non-retarded forces contribute 
to Tdis(rf) which consequently can be written as 



F dis (d) = - 



iff 



127rd 2 ' 



(4) 



where Ah is the Hamaker constant for the given layered 
system [l|, 0|. For symmetric systems, such as the ice- 
liquid-icc configuration, which occurs in the case of grain 
boundary melting, Ah is positive and hence dispersion 
forces are attractive. At the interface between the ice and 
the porous matrix Ah can assume both signs, depending 
on the material properties of the porous medium. 

Immobile surface charges with a given surface charge 
density q s on both interfaces which bound the quasi- 
liquid layer are screened by the counterions in the liquid. 
Making use of the Debye-Hiickel theory this leads to the 
contribution 



F e lec(d) 



< 

Kee 



(5) 



where eo is the vacuum permittivity, e the relative per- 
mittivity of liquid water, and 



' eepfc/jT 
e 2 N APl 



eeoksTd 
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(6) 



is the Debye length for monovalent ions, with the electron 
charge e, Avogadro's number Na, and Boltzmann's con- 
stant Ub- The Debye length describes the characteristic 
decay of the ion field adjacent to the charged surface due 
to the screening by counter ions. A repulsive force be- 
tween two charged surfaces originates in the restriction 
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of the entropy of the ions as the surfaces are brought 
closer. Whereas ions of charge opposite to those of the 
surface are attracted to it, they are repelled from each 
other, and the increased proximity induced by decreas- 
ing the film thickness increases the free energy. Note the 
relatively complex d-depcndcnce of F e \ ec (d) which is due 
to the fact that kT x is also a function of d. Thus, the 
richness of the role of impurities in the premelting of ice 
is a consequence of this and the efficient segregation of 
ions into the film; since nearly all of the impurities remain 
in the liquid, as the temperature changes the amplitude 
and range of the repulsive electrostatic interaction varies 
Q. In principle, the situation could be even more com- 
plex when for example q s = q s (d). Indeed, when d is 
small the Debyc-Hiickel expression is less accurate and 
stcric effects become more important. However, to intro- 
duce such complications here exceeds our understanding 
of the situation based on present experimental evidence. 

Using these expressions for 1(d) = Fdi S (d) +F e i ec (d) in 
Eq. ([3]) we obtain an implicit equation for the thickness d 
of the premelted layer as a function of the undercooling 
AT: 



PlQr 



AT 

T 



RgT m Ni 



A 



d 6ird 3 
where the constant c is 



-SL 



/N~3 

(7) 



e 2 N A 



y ee a k B T m 

This result is analyzed in detail in Section fill Al 



defines the surface normal. Thus, the surface energy, it- 
self a Gibbs free energy conjugate to area, is known in 
general to depend on the size of the system, the adsorp- 
tion of impurities and the crystallographic orientation |l| 
but the simplification of treating it as a constant is jus- 
tified both on the basis of experimental evidence for the 
first two dependencies and on the level of our analysis. 
The condition for thermodynamic equilibrium is again 



dG 



dG 

m 



(10) 



and hence a calculation analogous to that in the Sec- 
tion QTX] yields the result 



(ii) 



2jsi 

T m pi p s r 

or, for the general case of a surface with principal radii 
of curvature r\ and T2, 



AT pi _ . 



Pl 



Pa 



1 

'2 



(12) 



This result is the well known Gibbs- Thomson effect, i.e. 
the freezing point depression due to convexity of the ice- 
liquid interface, calculated here for the case where im- 
purities are present in the liquid water. The larger the 
undercooling AT, the smaller the radii of curvature of an 
ice-liquid interface in equilibrium. At the bulk freezing 
point, AT \ R g Tl Pl / q m pi, the radius of curvature di- 
(8) verges and a planar interface is the stable configuration. 



C. Premelting in a porous medium 



B. Curvature induced melting: The 
Gibbs-Thomson Effect 



In order to assess the effect of the ice-liquid interfacial 
curvature on the freezing temperature, we neglect the 
known anisotropy of the surface energy of ice [l[ and 
consider a spherical ice crystal of radius r surrounded by 
undercooled water. As opposed to the planar case, the 
interfacial area is no longer a constant as it was in the 
treatment of Section lll Al With this new complication we 
reconsider the Gibbs free energy G of the total system. 
The same holds for the numbers of moles N s , Ni, and Ni 
which are now taken for the total system and not per unit 
surface. For ease of notation we use the same symbols for 
these new quantities. All other quantities are as defined 
above. Whence, we have 

G(T,P,N s ,N h Ni) 

= p s (T, P)N S + p t (T, P)Ni + pi (T, P)Ni 



-RgTN, 



In — - 1 

N 



(9) 



where r is related to N s through N s — ^-r 3 p s with the 
constant molar density p s . We note here that we under- 
stand that 7 S ; = j s i(r,Ni,(fi), where <j> is the angle that 



We use the random close packing (RCP) of hard 
spheres as a model for the porous matrix. Apart from 
choosing the matrix material and hence calculating the 
properties which enter the determination of the Hamaker 
constant (see Section lll A [I . the only free parameter char- 
acterizing the porous matrix which contains the partially 
premelted ice is the radius R of the hard spheres. The 
properties of RCP hard spheres which are required for our 
calculation are the sphere packing fraction -q and the av- 
erage coordination number C. Recent estimates of these 
numbers are 77 = 0.634 and C = 6. It should be noted, 
however, that the basic concept of RCP is not well de- 
fined mathematically, so these numbers vary slightly ac- 
cording to the protocol and algorithm used to determine 
the exact configuration Q. 

At a given undercooling AT we are interested in the 
fraction of the interstitial volume which is occupied by 
quasi-liquid water fi which we will more compactly refer 
to as the liquid fraction. The first contribution we con- 
sider comes from interfacial premelting at the surface of 
the porous matrix. In the limit where the premelted layer 
is much thinner than the radius of the matrix particles, 
d <C R, this fraction is well approximated by 



/r at (AT, P . t ) 



^mat 



xd(AT, Pi ; AT, (13) 
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where a ma t is the surface area of the matrix particles 
per unit volume and d is obtained from Eq. ([7]) using 
the Hamakcr constant A^ at for the ice- water-matrix lay- 
ered system and the appropriate surface charge density 
q™ at . The sphere number density of the matrix is n s = 



r\ I '(47r R? 1 '3) and thus we obtain a ma t 
and 



4?ri?V = 3r]/R 



fr\AT, Pl ) = 



3y 

1-7/ 



xd(AT,pi;A^,g^)R- 1 . (14) 



The second contribution to the liquid fraction is due 
to grain boundary melting. Its importance depends sen- 
sitively on the grain size of the interstitial ice, which is a 
function of such factors as (1) the rate of cooling, (2) the 
undercooling AT, and (3) the time allowed for annealing. 
Even were we to attempt a complete computation of the 
coarsening rate of the ice grains in the matrix, there is 
no unique way to specify the initial value problem and 
hence to arrive at a concrete grain size as a function of 
the properties of the porous matrix. Thus, we introduce 
as an additional parameter the grain boundary density 
p g b, which is the surface area of grain boundaries per 
unit volume of the interstitial ice. We denote p g b = p g bR 
as the related dimcnsionless density. The corresponding 
liquid fraction follows immediately 

if laycr (AT, Pl ) = p gb x d(AT, Pl - <, qf)R- 1 , (15) 

where is the Hamaker constant for the ice-liquid-ice 
layered system and qf b is the surface charge density for 
the grain boundary. The thickness d of the quasi-liquid 
layer is again obtained from Eq. (0 . 

We now deal with curvature induced premelting. Im- 
portant contributions to the liquid water fraction come 
from regions where the ice-liquid interface has a large 
curvature. In our idealized matrix these are the points 
where two matrix spheres are in contact. Following CDF 
Q there is a liquid pocket forming at the contact point 
with a volume 2irr 2 R to lowest order in r, which is the 
radius of curvature of the ice-liquid interface (Fig. [T|). 
Taking into account that the second principal radius of 
curvature is much larger than r, we can use Eq. (|12|) to 
relate r to the undercooling AT: 



( \rp \ 1st ( AT p l 

r{AT,pi) = — q m — R g T m — 

Pa \ J-m Pl 



(16) 



Using the number density n s of the sphere matrix and 
the matrix coordination number C we can express the 
corresponding fraction of liquid water with respect to the 
total interstitial volume: 



jf™ (AT, Pi ) 



C/2xn s x2ir[r{AT, Pl )} 2 R 

i-n 

3 f] C \r(AT, Pi ) 

- X X — X — 

2 1-7] 2 [ R 



(17) 




FIG. 1: Sketch of the liquid pocket which forms between two 
matrix spheres of radius R. The amount of premelted liquid 
depends on the characteristic radius of curvature r of the ice- 
liquid interface as obtained from Eq. (|16|) . Depending on the 
matrix material there is either always a premelted liquid layer 
(repulsive dispersion forces, Ah < 0) or the premelted layer 
may collapse under certain conditions (attractive dispersion 
forces, Ah > 0) and the liquid-ice interface forms an angle 9 < 
tt/2 with the matrix surface. For parsimony of presentation 
we show the neighboring matrix spheres of different materials 
while the present theory assumes porous media made up of a 
single material. 



A further curvature induced contribution to the liquid 
fraction brought out by CDF comes from the lines where 
a grain boundary and a matrix sphere meet. 

In the limit that the matrix particles are small on the 
scale of the linear scale of the ice grains, but that R 
is much larger than r(AT, P i) the corresponding liquid 
fraction can be calculated using elementary geometry (see 
Appendix [A} 



/f bmat (AT, Pl ) 



= 3 | 2 - — I x x p gb x 



1-7? 



jAT, Pi ) 
R 



(18) 



Finally, curvature induced premelting gives rise to liq- 
uid veins along the lines where several grain boundaries 
meet. With the simplifying assumption that the ice 
grains arc arranged in a cubic packing, the total length 
of veins per interstitial volume, or the density of veins 
Pvoins, is obtained from p gb as /9 VO in S = Pg b / 3 - Given 
that the cross-sectional area of a vein in cubic packing is 



(4 



, we can calculate the corresponding contribu- 



with r(AT, pj) from Eq. ([T 



tion to the liquid fraction as 
/ ( veins (AT, pi) = Pvcins x (4-tt) x [r(AT, Pi )} 

= i(4-.)x&x(?*^] .(19) 

Therefore, the total liquid fraction is obtained as the 
sum of the five individual contributions which we have 
derived above; /, mat ,/f blayor j /contact^ /s bmat ? an d /™». 



D. Implicit equation for the total liquid fraction 

Thus far we have taken impurities into account in 
terms of a molar density pi or, in the case of the pre- 
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melted layer, a surface density of impurities iVj = pid. 
Laboratory probes using scattering or related approaches 
involve relatively small samples with good thermal con- 
trol and a well defined (fixed for a given run) dopant 
level The temperature changes in these systems are 
quasi-static so that the sample is at coexistence. In the 
natural environment there are of course temperature gra- 
dients. Therefore, for applicability of our work there, the 
typical broad temperature gradients found for example 
in natural soils imply that a local region within will be 
nearly isothermal and thus isosaline. Clearly, there will 
be many environmental circumstances where this may be 
violated, but because the first test of such a theory should 
be of high precision we focus mainly on the simplest ex- 
perimental setting and exploit the fact that the impu- 
rity segregation coefficient of ice is small Q so that the 
premelted liquid is enriched/diluted when the system is 
cooled/warmed. Hence it is prudent to introduce po, the 
molar density of impurities in the interstitial water when 
the system is above the bulk melting temperature. Thus, 
because of the validity of the assumption that impurities 
are contained almost exclusively in the liquid fraction of 
the interstitial water it follows that 



Po 

Pi = -r 
Ji 



(20) 



thereby properly accounting for the increase in the im- 
purity concentration in the premelted liquid as the liquid 
fraction decreases. Due to the fact that the regions of pre- 
melted water are all interconnected, pi is homogeneous 
over the entire sample. 

Having now considered the essential physical phenom- 
ena we can formulate an equation for the total liquid 
fraction /* ot in the porous medium as a function of the 
undercooling AT and the impurity concentration po in 
the completely melted sample. Combining all the contri- 
butions discussed in Section [II CI we obtain 



ft 



/r^AT, A,// Z tot ) + / ( s yer (AT, po/fn 



+ / co„tact (A7>o// to t) 
+ / vei„s (ATjp0// tot ) 



/f bmat (AT, P0 //r) 
(21) 



which is an implicit function for / ; tot = /; tot (AT, po) and 
in the general case has no analytical solution. Moreover, 
the underlying equation (JT]) for the thickness d of the 
premelted layer, which is required for /™ at and y ; gblayer ; 
is itself an implicit function of d. This calls for a care- 
ful stepwise analysis first of Eq. ([7|) and subsequently of 
Eq. (f2Tj) , which we perform in Sections MI Al and IIIIB1 



III. ANALYSIS AND APPLICATION 



A. Thickness of the premelted layer 



In Section III Al we derived Eq. Q for the thickness d 
of the premelted layer as a function of the undercooling 



AT, the impurity concentration JVj, the dispersion forces 
as embodied in the Hamaker constant Ah, and the sur- 
face density q s of immobile charges on the solid-liquid 
interfaces. This expression cannot be solved analytically 
in the general case. Hence, presently we analyze Eq. ([7]) 
analytically, to demonstrate that we can obtain known 
distinguished limits and facilitate intuition, and numer- 
ically, to provide specific ranges of results for particular 
physical systems of interest. We begin by rewriting this 
equation in terms of the reduced, or dimensionless, un- 
dercooling t = AT/T m as follows 



t = 



bN, 



a 

d? 



i - 



i 



/N~d. 



(22) 



where the constants b = R g T m /{piq m ), a = 
AH/(6irpiq m ), and q = q^ / (eeopiq m )- For attractive dis- 
persion forces which suppress premclting we have a > 0, 
and for repulsive dispersion forces which promote pre- 
melting we have a < 0. All other quantities in Eq. (|2"2")l 
are positive. 



1. Maximal undercooling and minimal layer thickness 

The simplest case occurs when both dispersion and 
electrostatic forces are absent, i.e., a = and q = 0. 
The thickness d of the liquid film is them simply given 
by 



d = bNd~ 



(23) 



representing the colligative effect of the dissolved impuri- 
ties which, due to the small segregation coefficient of ice, 
are strongly enriched in the liquid layer. The situation 
immediately becomes more complex if we take disper- 
sion forces into account. Let us first consider materials 
that have a > 0, which means that dispersion forces sup- 
press the formation of a premelted liquid film. Assuming 
q = for the moment, Eq. (f22j) possesses a solution as 
long as the undercooling t does not exceed a critical value 

< max = | — corresponding to a minimum film thick- 
ness dmin = y^jf-- If the temperature is reduced such 

that t > < m in the attractive dispersion forces dominate 
the colligative effect and the liquid film collapses to zero 
thickness. 

For the most general form of Eq. (|22[) . the maximal un- 
dercooling t max and corresponding minimal layer thick- 
ness dmin at which the quasi-liquid layer collapses cannot 
be expressed analytically. We obtain the global maxi- 
mum (dminjimax) by extremizing the right hand side of 
Eq. (|22|) with respect to d. Two limiting cases can be 
distinguished. For large impurity concentrations JVj the 
contribution due to surface charges becomes negligible 



and hence we recover the result 



3a 



2 (WV 4 ) 3 / 2 



/3a 



(24) 



for the case q = which we derived above. On the other 
hand, it can be shown that in the limiting case iVj — > 
the values d m i n and i max follow from the surface charge 
dominated limit (6 = and a = 0) which has 



dmm(iVi -> 0) 



2.62 



C (7V, -> 0) -> 0.0757 x g. 



(25) 

In Fig. [2] we plot the minimal thickness d m i n of the 
quasi-liquid layer corresponding to the maximal under- 
cooling i max as a function of the impurity concentration 
Ni for different values of the surface charge q s . We use a 
Hamaker constant Ajj = 3.0 x 10~ 23 J which is the re- 
sult obtained by Wilen et al. [3] for the ice/liquid/fused 
quartz layered system. The analytical results for Ni —> 
and Ni — > oo are characterized by the ^-dependencies 
of the asymptotic behavior. Note that the minimal thick- 
ness in the dispersion force dominated limit (q s = 0) is 
below 0.1 nm and therefore equal to the dimension of 
a single molecule. Hence, the notion of a quasi-liquid 
layer can no longer be applied and for the ice/liquid/fuscd 
quartz system layer collapse can only be observed exper- 
imentally at low impurity concentrations and it is always 
triggered by surface charges. This is seen directly from 
Eq. (|22]l where the contribution of the surface charge for 
sufficiently small d becomes negative, indicating attrac- 
tive forces, which suppress liquid layer formation. In or- 
der to observe a collapse of the quasi-liquid layer driven 
by dispersion forces, the latter have to be more strongly 
attractive. We have included in Fig. [2] the case of the 
ice-liquid-gold system which has a Hamaker constant of 
Ah = 1.573 x 10~ 21 J Q, and hence dispersion forces are 
indeed strong enough to yield a liquid-layer collapse at 
thicknesses of the order of those calculated for the sur- 
face charge dominated case. The contrast between the 
two scenarios will be studied further in the second part 
of this section (Figs.|4]and[5]). In Fig.|3]thc corresponding 
results for the maximal undercooling AT max = T m i max 
are plotted. As we expect from Eqs. (|2"4"]l and (f2"5"j). 
AT max increases with surface charge g s resulting in a 
g s -dependent offset for Ni — > 0. Moreover, for all q s , 
AT max increases with impurity concentration Ni . A large 
Hamaker constant, for example of the ice/liquid/gold sys- 
tem, results in a low AT max as shown by Eq. (f24| , where 

— 1/2 

AT max oc A H . There are interesting experimental con- 
sequences of Fig. |3l Obviously, measuring AT max for just 
one given impurity concentration Ni is not sufficient to 
determine the surface charge q s at the interface. This 
follows from the fact that any two curves intersect, and 
hence measurements for at least two different impurity 
concentrations are required. 



1 1 

ice/liquid/fused quartz system 

-I 




q s = 0.1\ 



ice/liquid/gold system, 
with q = 



\ q s = 0.5 q v = 0.75 




FIG. 2: In the case where dispersion forces at the interface 
are attractive, the quasi-liquid layer always collapses discon- 
tinuously from a minimal thickness d m i n , plotted here as a 
function of the impurity density Ni for different values of the 
surface charge q s , to zero when the temperature is lowered 
such that the undercooling exceeds AT max (plotted in Fig.[3|. 
The collapse is caused by the combined effect of immobile 
surface charges (dominating as Ni — > 0) and dispersion forces 
(dominating as Ni — >• oo). For the ice/liquid/fused quartz sys- 
tem the collapse caused be dispersion forces occurs at thick- 
nesses on the scale of the size of a single molecule where the 
notion of a quasi-liquid layer no longer applies. In order to ob- 
serve collapse induced by dispersion forces a larger Hamaker 
constant is required such as provided by the ice/liquid/gold 
system. 



2. Layer thickness for a given undercooling 



When q = 0, Eq. (|22|) can still be solved analytically 
but because the resulting expression is lengthy and not 
particularly instructive we present the solution only in 
the small t limit of moderate undercooling, 



b_Ni _^L_ 0(t3] 
t FN? 1 ' 



(26) 



This result clearly illustrates the fact that the presence 
of attractive dispersion forces (a > 0) acting across the 
system suppresses the formation of a quasi-liquid layer, 
whereas repulsive dispersion forces (a < 0) tend to in- 
crease the thickness of the premelted layer. 

When surfaces charges on the interfaces of the quasi- 
liquid layer with the substrate/ice are involved, Eq. (j2"2")) 
becomes transcendental and has no analytical solution. 
In order to get some insight in the effect of surface 
charges, we consider several special cases. First, let us as- 
sume that for moderate undercooling the system is dom- 
inated by surface charges. This scenario corresponds to 
setting a = 0, b = in Eq. (|2"2"j) and considering only 
the leading contributions as d approaches infinity. The 
solution is then 



[ln(f)] S 
c 2 N 



(27) 
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FIG. 3: The maximal undercooling AT max as a function of the 
impurity concentration Ni for the same parameters as used for 
the calculation of d m i n in Fig. [21 The zero slope asymptote 
for Ni —¥ is observed as expected from Eq. (|25[) . The large 
Hamaker constant of the ice-liquid-gold system leads to low 
values of AT max . 



This result reveals an interesting consequence of the pres- 
ence of charges at the interface. Under appropriate con- 
ditions the layer thickness d decreases when the impurity 
concentration increases. This observation is contrary to 
the naive expectation based on the colligative effect and 
originates in the fact that for larger N the electrostatic 
repulsion is screened for effectively [f| ■ 

For repulsive dispersion forces (a < 0) the solution of 
Eq. (|22|) can be obtained asymptotically in the limit of 
large undercooling t — > oo, whereas for attractive dis- 
persion forces (a > 0) the premelted layer collapses to 
zero thickness for t larger than some finite value. The 
repulsive case is 



d(t,N;a,q) 



a/3 



t l/3 



bN 
3i 



(-a)^q 



(28) 

and seems to contradict Eq. (f2"Tf for two reasons. First, 
the negative sign of the term proportional to q in Eq. (|28|) 
indicates that immobile surface charges impede the for- 
mation of the premelted layer, whereas in Eq. (|27p im- 
mobile surface charges alone are sufficient to stabilize a 
premelted layer. Second, while Eq. (|27fl implies that d de- 
creases when the impurity concentration Ni is increased, 
the opposite conclusion follows from Eq. ([2"8]l where, due 
to the term oc q, the inhibition due to surfaces charges is 
reduced when Ni increases. 

The apparent contradiction can be reconciled by the 
following observations. A direct consequence of Eq. (j2"2")l 
is that immobile surface charges do not contribute to the 
film thickness d when the latter equals the Debye length, 
namely do = c i\f. > i n which case the surface charge term 
vanishes. The film thickness d Q translates into an under- 
cooling to = bc 2 N 2 — ac^Nf . For an undercooling t < to 
surface charges enhance the formation of a quasi-liquid 



layer while for t > to surface charges suppress layer for- 
mation. Considering that Eq. (|2"7|) applies to the t — > 
case while Eq. (|28j) to the f ^ oo case, this explains the 
essential difference between Eqs. (f2~T)) and (j2"8")l in terms 
of whether they support or suppress the formation of a 
liquid film. 

In order to isolate the dependence of d on the impu- 
rity concentration Ni we consider Eq. (f2"2")l in the case 
where dispersion forces can be neglected (a = 0) and the 
immobile interfacial charge density is small (i.e. q — > 0). 
To lowest order d depends linearly on q, and hence we 
obtain 

d(t,N;a = 0,q->0) 



& f + ^-^)x e W^ + e%2 ).( 29 ) 



From this result the sign of dd/dN can be shown to equal 
the sign of 

S = a 2 {2-a) + e a /f3, (30) 

where a = cN l ^fbji and (3 = q/(bc 2 Nf). We find S > 
for all a > as long as /3 < /3 C = e 4 /32 ~ 1.71. There- 
fore if N > \J qj ((3 c bc 2 ) the layer thickness d increases 
with N for any given undercooling t. For lower impu- 
rity concentrations (j3 > /3 C ) there is a finite interval for 
a which contains a c = 4, where S < 0. Therefore, at 
sufficiently low impurity concentration there exists an in- 
terval of intermediate undercooling t where d decreases 
when Ni increases. This result implies that for large un- 
dercooling t the layer thickness d must always increase 
with Ni which was precisely a conclusion obtained from 
Eq. ([28]) . However, the result cannot be applied to the 
charge dominated case, Eq. (|2"T)) . because the correspond- 
ing limit b — > renders Eq. ((2l?)) trivial. It is possible 
though to obtain S < in the related limit b — >• 0, t — >• 0, 
if we assume that b/t approaches a large constant such 
that a > 2 holds. 

Attractive dispersion forces The non-monotonic de- 
pendence of d on Ni is illustrated in Fig. |4] where the full 
numerical solution of Eq. (f22j) is plotted for the case of 
the icc-liquid-fuscd quartz layered system. At low under- 
cooling AT the system displays characteristic colligative 
behavior as described by Eq. ((2"3"]) . The effect of immo- 
bile surface charges vanishes in the limit t — > as in- 
dicated by Eq. (f2U]) and the layer thickness d increases 
when the impurity density N% increases. The layer thick- 
ness for small AT is obtained by requiring that the im- 
purity concentration pi = Ni/d equals the impurity con- 
centration which would be necessary in a bulk system to 
see a freezing point depression equal to the undercooling. 
When AT is increased, surfaces charges start to increase 
the layer thickness relative to the system without sur- 
face charges (q s = 0). Interestingly, the resulting layer 
thickness for the system with the lower impurity density 
(N = 0.5 ^mol/m 2 ) becomes even larger than for the 
system with more impurities (N = 2.0 /imol/m 2 ). This 
illustrates the result from the above analysis displaying 
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FIG. 4: Thickness d of the quasi-liquid layer for the 
ice/liquid/fused quartz system (Ah = 3.0 x 10~ 23 J) as a 
function of the undercooling AT. For moderate undercool- 
ing the system is dominated by the colligative effect. In this 
regime the layer thickness increases with impurity density Ni . 
At intermediate undercooling the effect of immobile surface 
charges leads to a thickening of the quasi-liquid layer which 
is most effective for low impurity concentrations leading to a 
decrease of the layer thickness when the impurity density is 
increased. At large undercooling surface charge induced layer 
collapse is observed, while layer collapse due to dispersion 
forces alone occurs at layer thicknesses below 0.1 nm. 



the existence of an interval of intermediate undercooling 
where the layer thickness decreases when the impurity 
density increases. Recall that the reason for this is that 
at low impurity concentrations the electrostatic force is 
less effectively screened and hence the repulsion due to 
surface charges has a stronger effect, leading to a larger 
layer thickness than for higher impurity concentrations. 
When the system is undercooled further the collapse of 
the quasi-liquid layer can be observed in the figure. Here, 
because of the weakness of dispersion forces, the collapse 
is a consequence of the electrostatic force becoming at- 
tractive at low layer thicknesses. Collapse for the cases 
without surface charges cannot be observed in Fig. @] as 
it occurs at thicknesses below 0.1 nm (c.f., Fig. [2]). In 
order to see layer collapse at larger layer thicknesses a 
larger Hamaker constant is required, as provided by the 
ice-liquid-gold system. To see this, we plot solutions 
of Eq. (|22|) for this system in Fig. [5] for the same pa- 
rameters as considered in Fig. 2] for the ice-liquid-fused 
quartz system. Layer collapse is now observed at thick- 
nesses above 0.1 nm even for the configurations with- 
out surface charges. For the configurations with surface 
charges we observe that for Ni = 0.5 /miol/m 2 layer col- 
lapse is not altered by the stronger dispersion forces while 
for N = 2.0 /umol/m 2 the stronger attractive dispersion 
forces lead to a layer collapse at a larger layer thickness 
compared to the corresponding configuration considered 
in Fig. 21 Thus, in this case stronger dispersion forces 
drive layer collapse. 
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FIG. 5: Thickness of the quasi-liquid layer for the 
ice/liquid/gold system (Ah = 1.573 x 10~ 21 J). In contrast 
to the ice/liquid/fused quartz system, Fig. [4j layer collapse 
occurs at thicknesses above 0.1 nm even for the cases with- 
out surface charges. This is due to the stronger attractive 
dispersion forces in the ice/liquid/gold system. Layer col- 
lapse for the cases with surface charges is altered only for 
the higher impurity density Ni = 2.0 /imol/m 2 relative to the 
ice/liquid/fused quartz system, because at d ~ 1 nm disper- 
sion forces are weak compared to the electrostatic force. 



Repulsive dispersion forces The result of the analysis 
following Eq. (f29|) can be tested on the ice/liquid/silicon 
system which has a Hamaker constant of Ah = —1.66 x 
10~ 21 J [H, the negative sign implying repulsive disper- 
sion forces. As shown in Fig. [6] the system displays 
the expected colligative behavior for low undercooling 
and an enhanced layer thickness for intermediate un- 
dercooling where the effect of surface charges enhances 
the colligative effect, which alone would yield a straight 
line, d cx N^~ . For intermediate undercooling the fig- 
ure reveals a transition from a regime of low impurity 
density where the layer thickness decreases when N 
increases, to a regime of high impurity density where 
the colligative effect dominates and accordingly an in- 
crease of d is observed when N increases. From the fig- 
ure the locus of the transition can be estimated to be 
AT C « 10 - 20 K, d c « 1 nm, and N i>c w 2.5 /imol/m 2 . 
This is in good agreement with the analytical result ob- 
tained from Eq. AT C = 9.6 K, d c = 1.25 nm, and 
Ni c = 3.97 ^mol/m 2 . It should be noted that the an- 
alytical result is exact only in the limit q s — > but the 
complete solution of Eq. (f22|) shows that it applies quali- 
tatively and semi-quantitativcly for surface charge densi- 
ties as high as q s = 0.5 C/m 2 . For the ice/liquid/silicon 
system no layer collapse is observed. This results from 
the fact that dispersion forces, which dominate at low 
layer thicknesses, are repulsive and stabilize the quasi- 
liquid layer. However, at impurity densities Ni smaller 
than those considered in Fig. [6j the Debye-Hiickel ex- 
pression for F c \ cc generates an attractive force in the 
ice/liquid/silicon system with q s = 0.5 C/m 2 at suffi- 



9 



10 



0.1 



N. = 1.25 nmol/m 

N. 



i 2.5 |Xmol/m 

2 

■ 5.0 umol/m 



q s = 0.5 C/m 

d =1.25 nm 

c 

N. = 3.97umol/m 2 



AT = 9.6 K \ s s. 



0.1 



10 



100 



AT [K] 



FIG. 6: Thickness d of the quasi-liquid layer for the 
ice/liquid/silicon system (Ah = — 1-66 x 10 -21 J) as a func- 
tion of the undercooling AT. As dispersion forces in this sys- 
tem are repulsive they stabilize thin quasi-liquid layers and 
hence no layer collapse is observed. The figure illustrates the 
analytical result obtained from Eq. ([29} which states that 
for intermediate undercooling AT and Ni below a transition 
value the layer thickness decreases when the impurity concen- 
tration increases. For large Ni the colligative effect dominates 
and d increases when Ni is increased. The transition value 
Ni, c as well as the locus AT C and d c where the interval with 
dd/dNi < appears as derived from the figure are in good 
agreement with the analytical result (values given in the fig- 
ure). Note the sigmoidal transition to small but finite film 
thickness for the lowest value of Ni. 



ciently small layer thicknesses which results in a discon- 
tinuous drop in the layer thickness at a certain under- 
cooling. The layer, however, does not collapse to zero 
thickness but rather is stabilized at a smaller thickness 
due to the repulsive dispersion forces which dominate for 
very thin quasi-liquid layers. In this article we do not 
investigate this phenomenon of layer discontinuity which 
occurs in a regime where both the d-dependence of q s , 
the peculiarities of the Debye-Hiickel result, and other 
effects, have to be accounted for. 



B. Volume fraction of premelted water 

Having analyzed prcmelting in planar layered systems, 
we now return to the implicit equation (Eq. I21[) for the 
volume fraction of premelted liquid at a given undercool- 
ing t that we introduced in Section III Dl In addition to 
being implicit for /; the equation uses the thickness of the 
premelted layer, which itself is obtained from an equation 
which does not have an analytical solution (Eq. |2"2")) . Ac- 
cordingly, we do not possess an analytical expression for 
/; in the general case. 

Low surface charge However, where immobile surface 
charges can be neglected (q s — 0), Eq. ([2"Tj) becomes more 
tractable. In the case where the matrix material has a 



negative Hamaker constant (A^ at < 0), i.e., dispersion 
forces are repulsive and support the formation of a pre- 
melted layer, Eq. ([2"Tj) becomes 
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with £ = j s i/p s q m and where c = R g T m p /q m pi is, up 
to a constant factor, the density of impurities po in the 
fully melted system. 

Obviously, Eq. (|3"Tj) requires that t > Co/// must hold 
and thus, because < fi < 1, when t < cq Eq. (|3"1| 
is not applicable. The case t < cq corresponds to an 
undercooling which is smaller than the depression of the 
bulk freezing point due to the colligative effect originating 
from the presence of the impurities. Therefore, we must 
assume a liquid fraction /; = 1 when t < cq. 

While Eq. (|3"Tj) does not posses a general analytical 
solution, it can be solved in the limit of low impurity 
concentration Cq to yield 
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-0(c§), (34) 

(35) 
(36) 



Note that the term oc Co is sufficient to obtain the cor- 
rect leading contribution due to impurities in both lim- 
iting cases t — !- and t — >• oo. The factor oc c\ becomes 
a constant in the limit t — > where vanishes. The 
coefficient belonging to Cq vanishes as i~ 7 / 3 in the limit 
t — > oo and is thus subdominant to J| . 

In the case where dispersion forces are attractive 
(A 1 ™' > 0) the formation of a premelted layer requires 
an impurity concentration which is higher than the one 
which would be required to prevent the bulk liquid from 
freezing. Hence, in this case all the liquid water in the 
matrix pores must originate from liquid veins which form 
due to curvature induced premelting. In contrast to the 
premelted layer with Af^ > these veins have formed 
when t > Co//;, a temperature for which the premelted 
layers are still collapsed and do not contribute to the 
liquid fraction. Thus, premelted layers do not actually 
form because at higher temperatures, beginning with the 
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veins, the sample completely liquifies. Therefore, the cor- 
responding equation for // is 



fl 



D 



(37) 



which has a simple analytical solution. Again, we exam- 
ine the leading order contributions in the limiting cases 
t — > and t — > oo to find 




where the latter relation is an extension of Eq. which 
is undefined in the limit t — > oo and A — > 0. 

Non-zero surface charge In the general case with im- 
mobile surface charges (q s ^ 0) we solve Eq. (|2~Tj) numer- 
ically. In Fig. [7] we show results for a matrix of fused 
quartz with a grain size of R = 500 pm, which may serve 
as a model for frozen soils. The surface charge density 
is chosen to be q s = 0.5 C/m 2 which has the correct 
order of magnitude for mineral/ water interfaces at high 
pH [13] ■ The ice which fills the voids of the porous ma- 
trix is taken to be polycrystalline with a grain bound- 
ary density p g b = 1 pm -1 . A comparison of Eqs. (| 14[) 
and (|15|) shows that for the given pgb the surface area 
provided by grain boundaries is by a factor 100 larger 
than the surface area of the porous matrix. Therefore, 
we expect the impurity concentration in the sample to 
be very low when grain boundaries are melted and, ac- 
cording to Eq. (|24f , the maximal undercooling AT max for 
which quasi-liquid intcrfacial water at the grain bound- 
aries is stable should be low. An upper limit for AT max 
can be calculated by assuming that all the impurities in 
the sample are located at the premelted grain boundaries. 
For the values of po from Fig. [7] we use Ni = p / p g b to 
obtain a range jVj = 0.01 — 0.1 /xmol/m 2 and, with the 
aid of Eq. (|24|) this corresponds to a range in AT max of 
0.01 — 0.26 K. The range of AT max uses the Hamakcr 
constant Ah = 3.3 x 10 -22 J computed from complete 
dispersion force theory Q in the limit of thin layers with 
the data for the dielectric functions from Elbaum and 
Schick [111 ], Including a typical surface charge density 
of q s = 0.01 C/m 2 for the grain boundary [12j lowers 
these values to ATmax = 0.01 — 0.20 K. Hence, the grain 
boundaries are always collapsed to zero film thickness 
for the range of AT considered in the figure. A promi- 
nent feature of fi = ff ot with po = 0.01 mol/m 3 in the 
fused quartz matrix is the discontinuity at AT = 22.6 K 
where /; drops by almost an order of magnitude. This 
is the signature of the collapse of the premelted layer at 
the ice/fused quartz interface described in detail in Sec- 
tion QlT2] In the strongly undercooled, layer collapsed 
system, the only contribution to fi comes from curva- 
ture induced melting; /f urv = / ; c ° ntact + /s bmat + /veins 
(cf. Eq. |2~T|) , which discontinuously increases at the un- 
dercooling where the premelted layer collapses. This is 
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FIG. 7: Liquid fraction /; for water in a porous matrix of 
fused quartz with a grain size R — 500 pm as a function of 
the undercooling AT. The result for the low impurity density 
po = 0.01 mol/m 3 (with respect to the interstitial volume) 
shows a discontinuity where the quasi-liquid layer between 
the ice and the matrix collapses. This is accompanied by 
an increase of the impurity concentration in the remaining 
(curvature melting induced) liquid which creates a discon- 
tinuous increase of /; culv . For the higher impurity density 
po = 0.1 mol/m 3 the discontinuity is not visible in the figure 
because it occurs at large undercooling (AT « 110 K). 



because, due to the efficient rejection of impurities by 
the ice lattice, the impurity concentration in the remain- 
ing liquid increases discontinuously when the liquid layers 
collapse. At low undercooling curvature melting becomes 
increasingly important exceeding the contribution from 
interfacial premelting /' aycr = / ; mat . However, at larger 
undercooling, close to the location of layer collapse, cur- 
vature induced melting yields only a minute contribution 
to fi. The importance of curvature induced melting at 
low undercooling is reflected in Eq. (|34[) where the cur- 
vature term B/t 2 dominates for t — >• 0. In the opposite 
limit of large undercooling the analytical result Eq. ([55)1 
applies. The leading term co/t originates from the col- 
ligative effect and curvature induced melting enters only 
through the subdominant term. Hence, for large under- 
cooling the principal contribution to the liquid fraction 
is the colligative effect. However, liquid water is located 
in regions of high curvature; contact of matrix spheres, 
contact lines of grain boundaries and matrix spheres and 
liquid veins. Figure [7] also shows // in the same system 
with a higher impurity concentration (po = 0.1 mol/m 3 ). 
Here the collapse of the liquid layer is moved to a lower 
temperature (AT w 110 K) so that no discontinuity is 
observed in the range of undercooling considered in the 
figure. 

In Fig. [5] many of the same conditions as in Fig. [7] 
are considered (R = 500 pm, q s = 0.5 C/m 2 , p = 
0.1 mol/m 3 ) the principal difference being that the pore 
ice is single-crystalline (p g b = 0). We find that ff ot in 
the fused quartz matrix is qualitatively the same as for 



11 



the polycrystalline ice (see Fig. [7J) in that it does not 
show liquid layer collapse on the scale of the plot for 
po = 0.1 mol/m 3 . When replacing the matrix with gold, 
which has stronger attractive dispersion forces at the ice- 
gold interface, we see a collapse of the premelted layer at 
AT «7K. This is accompanied by, as in the case of layer 
collapse for the fused quartz matrix at po = 0.01 mol/m 3 
(Fig. [7]), a discontinuous increase of the liquid fraction 
/j curv . As for the fused quartz matrix, for sufficiently low 
undercooling that the premelted layer has collapsed, the 
only contribution to / ; tot comes from curvature induced 
melting; ff ot = /f urv . While, on the scale of the plot 
ff ot does not display a discontinuity at the point of layer 
collapse, there is a discontinuity A/; which can be seen 
by expanding the full result in the limit of large impurity 
concentration cq 

where t co \\ is the undercooling at layer collapse and /; Iayer 
the liquid fraction from interfacial melting at a temper- 
ature slightly above collapse t = t co i\ — St, with St > 0. 
Clearly, A/; — > for cq ^ oo. 

An important result seen in Fig. [5] is that although 
the Hamaker constants for gold and fused quartz differ 
by almost two orders of magnitude, the results for /; are 
almost identical over a large range of AT; the deviation 
increases from 0.02 % to 0.4 % at layer collapse and is 
6 % for AT = 60 K. This is because dispersion forces are 
only relevant for thin premelted layers and hence at tem- 
peratures slightly above layer collapse. However, when 
layer collapse occurs, according to Eq. (|4TJ)) the high im- 
purity concentration, combined with single-crystallinity 
(i.e. small B). insures that the Debye length and the 
discontinuity of // are small. Therefore, the deviation 
between the results for the two different matrix materi- 
als remains small. The conclusions drawn from Fig. [5] 
are important to note for experimental determination of 
// in different materials because they imply that even 
though the underlying microscopic behavior of two sam- 
ples may fundamentally differ (interfacial premelting vs. 
collapsed premelted layer) this may not be resolved in the 
integrated/homogenized quantity //. Spatially resolved 
measurements may be required to detect characteristic 
transitions such as the collapse of interfacially premelted 
layers. 

Effect of polycristallinity Finally, we study the in- 
fluence of the degree of polycrystallinity of the inter- 
stitial ice for the case where the porous matrix con- 
sists of silicon, which leads to dispersion forces that are 
repulsive (Ah < 0) and prevent the interfacially pre- 
melted layer from collapsing. As in the previous exam- 
ples the grain size of the matrix is R = 500 /xm and as in 
Fig. [7] the grain boundaries are collapsed. Thus, Fig. [S] 
shows the curves for /; = /' ot , the interfacial melting 

(// ayer ) an d curvature induced melting (/f urv ) contribu- 
tions for single-crystalline (p g b = 0) and polycrystalline 
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FIG. 8: Liquid fraction /( for the gold matrix compared to 
the fused quartz matrix as a function of the undercooling AT 
for po = 0.1 mol/m 3 . Importantly, although the Hamaker 
constants for the two materials differ by almost two orders 
of magnitude, which results in a layer collapse for the gold 
matrix at AT ss 7 K, the curves for the total liquid fraction 
are almost identical. Note that the curves for fused quartz 
have been shifted for clarity. In particular, the discontinu- 
ity of / ; tot at the location of layer collapse is too small to be 
resolved in the figure, which results from the high impurity 
concentration, see Eq. (|40|l . The strength of dispersion forces 
is important only for thin premelted layers, i.e. in the imme- 
diate vicinity of layer collapse. The discontinuity for fused 
quartz is not visible in the figure because it occurs at large 
undercooling (AT ^ 103 K). 



ice (p g b = 1 /im _1 ). As expected, the polycrystallinity 
leads to a higher liquid fraction In particular, for 
small undercooling AT, where curvature induced melting 
is enhanced and hence liquid veins greatly contribute, fi 
consists almost entirely of J ; curv . The values of / ; curv for 
the single crystal ice, for which contributions derive only 
from the contact points between matrix spheres, are al- 
ways by at least an order of magnitude lower than for the 
polycrystalline ice case. This is not as obvious as it may 
first appear. Due to the higher liquid fraction in polycrys- 
talline ice the impurity concentration pi is smaller than 
for the single crystal case. Thus, polycrystalline ice has a 
smaller ice- water radius of curvature (see Eq.[l6]). Conse- 
quently, the volume of liquid contained for example at the 
contact point between two matrix spheres, is smaller for 
the polycrystalline case than for the single crystal case. 
It is the grain boundary density p g b, which causes B to 
be larger for the polycrystalline ice, that insures J ; curv 
will be also larger, as indicated by Eq. In addition 

to the cases studied in Section fill A[ another interesting 
manifestation of the non-monotonic dependence of the 
premelted layer thickness on the impurity concentration 
for the systems with interfacial charges can be seen in 
Fig. [9] Although the impurity concentration pi and the 
impurity density iVj are always larger in the single-crystal 
case, the liquid fractions // ayor , which from Eq. ([T4"]) are 
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FIG. 9: Liquid fraction fi for a matrix consisting of silicon 
spheres with radius R = 500 /im as a function of the under- 
cooling AT. Two different cases concerning the ice, which is 
contained in the matrix, are considered: polycrystalline ice 
(p g b = 1 Mm -1 ) and single-crystalline ice (p g b = 0). The 
polycrystalline sample always has a higher liquid fraction, 
the difference with the single-crystalline sample being most 
pronounced at small undercooling. In this regime /; in the 
polycrystalline sample is almost entirely made up of / cmv , 
which mainly contains contributions from liquid veins in the 
ice. An interesting consequence of the non-monotonic depen- 
dence of the premelted layer thickness on the impurity den- 
sity Ni can be observed: although Ni is always larger in the 
single-crystalline sample than in the polycrystalline sample, 
the curves for /, laycl (which are proportional to the layer thick- 
nesses) intersect. 



proportional to the corresponding layer thicknesses, in- 
tersect at an undercooling AT int = 1.38 K. Thus for 
AT < AT; nt the larger impurity density of a single crys- 
tal sample leads to a smaller layer thickness than the 
polycrystalline case. It is only for AT > ATj nt that the 
usual increase of the layer thickness with increasing im- 
purity concentration is observed. This complex behavior 
provides further evidence for the observation (Fig. 
that measuring the integrated/homogenized quantity // 
alone docs not yield satisfactory information about the 
microscopic structure and behavior of the sample, in par- 
ticular the rich and interesting temperature dependence 
of the thickness of the premelted layer in a frozen matrix 
containing ice with different degrees of polycrystallinity. 
By parity of reasoning, inclusion of just one or two mi- 
croscopic phenomena in a theory that is intended to pre- 
dict the effective medium behavior of a porous medium, 
will provide irrelevant and misleading predictions unless 
the proper asymptotic limits of a complete theory are 
assured. 



IV. SUMMARY AND CONCLUSION 

In this work we construct a theory for the premclting of 
ice in a porous matrix (Scction|TI|, which is modeled as a 
random close packing of monodispcrse spheres. The the- 
ory includes premelting at the matrix-ice interface, grain 
boundary premelting, as well as curvature and impurity 
induced premelting. The latter includes premelting in 
liquid veins, at contact lines of the matrix with a grain 
boundary, and at the contact points between two matrix 
spheres. The matrix is characterized by the sphere ra- 
dius R and the matrix material, which enters the theory 
via the Hamaker constant Ah for the ice/liquid/matrix 
layered system. In addition to dispersion forces at the 
ice-matrix interface, which, depending on the matrix ma- 
terial, either stabilize or destabilize the premelted quasi- 
liquid layer, we include the effect of electrostatic forces 
at the interface due to adsorbed surface charges at the 
interface along the lines of our previous work Ref. 0. The 
theory includes the effect of charged impurities in the 
premelted water both on the thickness of the premelted 
layers (grain boundaries, ice-matrix interface) as well as 
on the curvature induced melting. The main idea of the 
present approach is to make use of the small segrega- 
tion coefficient of ice which leads to an enrichment of 
impurities in the premelted liquid (making the impurity 
concentration a function of the fraction of premelted liq- 
uid) and of the fact that regions of premelted liquid in a 
sample are all interconnected and therefore the impurity 
concentration is homogeneous. This approach allowed 
us to derive an implicit equation for the liquid fraction 
fi as a function of the undercooling AT, Eq. (|21[) . which 
we analyze for general parameters by considering limiting 
cases where analytical solutions are possible (Section lIII[) . 
To understand each contribution we began by perform- 
ing a careful study of intcrfacial melting, analyzing an 
implicit equation for the thickness d of the premelted 
layer as a function of AT, Eq. ([7]), for general parame- 
ters. In particular, we addressed the discontinuous col- 
lapse of the premelted layer under certain conditions and 
the non-monotonic dependence of d on the impurity con- 
centration (cf. Ref. [H). However, in the general case both 
Eq. (J?]) and Eq. (|2~T|) , the latter of which makes use of the 
former, have to be solved simultaneously through a nu- 
merical scheme. In so doing we chose physically relevant 
parameters of interest to a wide community of scientists. 
For the porous matrix we use R = 500 fim and considered 
two materials (fused quartz and gold) that yield attrac- 
tive dispersion forces at the ice/matrix interface, thereby 
inhibiting the formation of a premelted quasi-liquid layer, 
and one material (silicon) which enhances the premelted 
layer. The surface charge at the interface is taken to 
be q s = 0.5 C/m 2 . The impurity concentrations lie in 
the range of 0.01 — 0.1 mol/m 3 defined with respect to 
the interstitial volume, and we consider both single- and 
polycrystalline ice. We calculate /; for AT in the range 
of 0.5 — 60 K and obtain liquid fractions which range be- 
tween 10~ 6 and 10~ 3 of the interstitial volume. 
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For the parameters under consideration the grain 
boundaries are always of zero thickness (i.e., collapsed), 
and hence only the collapse of the quasi-liquid layer at 
the ice-matrix interface is observed in fi (see Fig [7] for 
the case of the fused quartz matrix). When the tem- 
perature of the sample is decreased below the point of 
layer collapse, only curvature induced prcmclting con- 
tributes to While dispersion forces are too weak in 
the fused quartz system to significantly influence the re- 
sults, in the case of the gold matrix dispersion forces are 
sufficiently strong to drive layer collapse, as opposed to 
the surface charge triggered layer collapse in the fused 
quartz system (see Fig. [8] for a comparison of the two 
cases). An important conclusion from Fig. [8] is that for 
high impurity concentration and/or low degree of poly- 
crystallinity the characteristic discontinuity of fi upon 
layer collapse is minute and most likely not resolved 
in experiments. Thus, transitions that would be dra- 
matic and easily detected at planar interfaces will be 
much more challenging to extract from measurements 
of the integrated/homogenized quantity //. Finally, we 
study the influence of the degree of polycrystallinity on 
fi for the case of the silicon matrix where, due to the 
repulsive dispersion forces, layer collapse does not occur 
(Fig. [5]). The results reveal an interesting consequence 
of the non-monotonic dependence of d on the impurity 
concentration. Although the actual impurity concentra- 
tion is higher in the single-crystalline sample than in the 
polycrystalline sample at all temperatures, the thickness 
of the premelted layer can be both larger and smaller in 
either case, depending on the temperature range. This 
peculiarity is obscured by the curvature induced contri- 
bution to /; providing further evidence that additional 
measurements are required // to capture the physics con- 
trolling the sample microstructure. 

We view our approach as having provided a thorough 
theory of premelting in a porous medium which can be 
applied to a plethora of experimental situations and ma- 
terial systems. For the sake of brevity and clarity, we 
presented examples of specific materials systems; gold, 
fused quartz, and silicon. However, we chose these ma- 
terials because their dielectric properties, as reflected in 
their Hamaker constants, span all the range of materials 
as diverse as sapphire, polystyrene, and polyvinylchlo- 
ride [1] . We did not explore the effect of varying the ma- 
trix sphere radius R, its roughness or its surface charge 
density q s . The value R = 500 /irn was chosen with the 
example of frozen soils in mind, while q s = 0.5 C/m 2 pro- 
vides a realistic order of magnitude as inferred from ex- 
perimental data for mineral-water interfaces at large pH 
[l0(. However, not much is known experimentally about 
q s in the ice/liquid/matrix systems and the simplifying 
assumptions of symmetry (same charge at the ice-liquid 
and liquid-matrix interfaces) as well as independence of 
q s on the impurity concentration and layer thickness are 
not yet established experimentally on planar systems. 
Nevertheless, given the importance of the impurity con- 
centration dependent screening length for regulating the 




FIG. 10: a) Geometry of a matrix sphere intersecting with 
a grain boundary which is located at distance z from the 
sphere center, b) Zoom on the curvature melting induced 
liquid "vein". As r <C R the matrix surface can be considered 
as flat in the vicinity of the liquid vein. 



strength and range of the electrostatic force, the precise 
value of q s may be considered secondary. 

An important question for future research concerns the 
dynamical evolution of the porous matrix system in re- 
sponse to changes in the control parameters, in particu- 
lar the sample undercooling AT. In order to guarantee a 
homogeneous impurity concentration throughout a sam- 
ple, AT must be changed on a time scale that is slow 
compared to the typical impurity diffusion time. An 
even more ambitious question concerns impurity diffu- 
sion during layer collapse. The situation where the pre- 
melted quasi-liquid layer between the ice and the matrix 
collapses, which forces the impurities to migrate to re- 
maining liquid regions with curvature induced melting 
constitutes a challenging task for further investigations. 
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Appendix A: Melting volume at a grain 
boundary—matrix contact line 

We consider a matrix sphere that is embedded at the 
interface between two ice grains. The grains are assumed 
to have a planar grain boundary and to be much larger 
than the matrix sphere. Let the distance between the 
sphere center and the grain boundary be z < R. Curva- 
ture induced melting leads to a circular "vein" along the 
line where the grain boundary and the sphere intersect 
fFig. UOb ,). Because r is assumed to be much smaller than 
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R we can consider the matrix surface as being planar in 
the vicinity of the vein as shown in Fig. [TUb . The angle 
between the grain boundary and the matrix sphere is 



(Al) 



a = arccos — 
\R 



Using elementary geometry we obtain the surface areas 
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Accordingly, the volume of premelted water for the con- 
figuration in Fig. [TU] is 2irhA to t . Taking all possible dis- 
tances of the matrix sphere to the grain boundary into 
account yields an average volume v of premelted liquid 
as 
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Now consider a random close packing with packing 
fraction r\ (the porous medium), which contains ice with 
two different crystallographic orientations being sepa- 
rated by a planar grain boundary. The volume fraction r/ 
is also the surface fraction of the grain boundary which 
is cut out in the form of circles with radii between and 
R due to the presence of the matrix spheres. Given that 
the average surface area per intersecting sphere is 



4 I" dzir(R 2 - z 2 ) = IttR 2 
R Jo 3 



(A6) 



a sample of volume V s = A s x s with a grain boundary 
perpendicular to the x-axis contains r\A s ja intersection 
circles of matrix spheres with the grain boundary, if A s 
is sufficiently large. If the sample contains the total grain 
boundary surface p'^Vs this corresponds to a net surface 
A s = p' h V s /(l — if) and hence a density of intersection 
circles 



VPgb = VPgb 
(1 — r/)d a 



(A7) 



where p g h is the grain boundary density with respect to 
the interstitial volume. The volume of premelted liquid 
along the intersection circles (per unit volume of the sam- 
ple) is vn c . Or, relative to the interstitial volume: 
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